On determining conditions and suitable locations for fish survival by using the solution of the two coupled pollution and aeration equations

The coupled equations of pollution and aeration for flow in a river were studied under generalized assumptions in terms of parameter dependency on space and time, as well as general boundary constraints. An analytical solution was obtained in the steady-state case. Also, the system was solved in its unsteady state numerically in a dimensionless form using the finite difference scheme. The effect of different parameters controlling the flow (such as the velocity, Peclet number, injected pollutants, and so on…) was studied. Investigations indicate that the special cases of the proposed model (i.e., uniform distribution of pollutant and Dissolved Oxygen concentrations, and zero injected pollutants along the river) give results that agree with the previous studies. This simple model helps in understanding the behavior of the pollution-aeration process and its relation to the injected pollution along a river and its effect on fish survival. A simple procedure was discussed in this study to help in regulating farming, industrial, and urban practices and impose restrictions if necessary. This study determines with accuracy the intervals of the river at which fish can survive at a given time, as well as the maximum amount of pollutants allowed to be injected along the river for fish survival.

Water pollution is one of the most important problems that directly affect human lives, fish survival, and environmental health in general. Hence, various attempts are made to search for affordable, efficient, and environmentally friendly water systems. Firstly, solar energy-based technologies 1 are proposed by researchers 2 as a cost-effective energy source in solar water pumps. Secondly, methods for heavy metal removal are being discussed 3 . Thirdly, biological wastes that can be treated by dissolved oxygen (DO) are modeled 4 , which is the focus of the current study. Predicting pollutant concentration in advance will reduce the negative effects drastically 5 . Tons of fish die every year as a consequence of water pollution 6 . Thus, mathematical models describing solute transport existed as early as the 1920s. Generally, the two main solute transport models are those based on the advection-diffusion equation (ADE), with different assumptions concerning initial and boundary conditions, the dependence of the parameters on space and time, and their generalizations in two and three dimensions. The other is the aeration model, which consists of two coupled non-linear partial differential equations (PDEs) with their respective conditions. The coupling occurs because the oxygen reacts with the pollutants that result from industrial and agricultural areas, producing harmless compounds. The aeration model is a powerful tool to help impose restrictions on urban and farming practices. The importance of this model is in enabling simulation scenarios to be tested for fish survival, which usually requires that the DO concentration be taken to be above 30% of the saturated DO concentration 7 . The pollution model has many analytical solutions for different cases. For instance, Shukla 8 discussed the unsteady transport dispersion of nonconservative pollutants with timedependent periodic waste discharge concentration. Jaiswal et al. 9 obtained analytical solutions for the temporally and spatially dependent solute dispersion of pulse-type input concentration in one-dimensional semi-infinite media. Manitcharoen and Pimpunchat 10 acquired analytical and numerical solutions of pollution concentration with uniformly and exponentially increasing forms of sources. Yadav and Kumar 11 solved the space and time-dependent ADE and obtained an analytical solution of two-dimensional conservative solute transport in a

The governing equations
The coupled non-linear PDEs modelling pollutant and DO concentrations in a river can be written as 4,10,11 : where R is the retardation factor, A(km 2 ) is the cross-section area of the river which is taken to be constant, is the diffusion coefficient of the pollutant concentration, D X (x, t) (km 2 day −1 ) is the diffusion coefficient of DO, u(x, t) (km day −1 ) is the river's velocity, K P (x, t)(day −1 ) is the degradation rate coefficient for pollution at a constant temperature, K X (x, t)(day −1 ) is the degradation rate coefficient for DO at a constant temperature, k(kg km −3 ) is the half-saturated oxygen demand concentration for pollution decay, q(x, t) (kg km −1 day −1 ) is the rate of the injected pollutant discharge, (km −1 ) is an arbitrary constant of pollution source terms, α(x, t) (km 2 day −1 ) is the mass transfer of oxygen from air to water and S(kg km −3 ) is the saturated oxygen concentration.
The flow domain parameters are considered to be spatially and temporally dependent so we have: where a (km −1 ) is a heterogeneity parameter and b is an arbitrary dimensionless constant 11 , R 0 is the initial retardation factor, D P 0 is the initial pollution diffusion coefficient, D X 0 is the initial DO diffusion coefficient, u 0 is the initial velocity component, q 0 is the initial rate of the injected pollutant discharge, K 1 is the initial degradation rate coefficient for pollution at a constant temperature, K 2 is the initial degradation rate coefficient for DO at a constant temperature, α 0 is the initial oxygen mass transfer, m (day -1 ) is a constant. We choose g(mt) such that for t = 0 or m = 0 we have g(mt) = 1 which ensures that the nature of the initial condition doesn't change in the new time domain T (day) which is defined by 24 : Equations (3 and 4) transform Eqs. (1 and 2) into: www.nature.com/scientificreports/ where U P = u 0 + aD p 0 , U X = u 0 + aD X 0 , q 1 = q 0 A ; α 1 = α 0 A . Assuming the initial solute concentration to be an exponentially decreasing function of the space variable 9,13 , a diurnal flow variation from a wastewater treatment plant with a constant discharge is modelled approximately by a sinusoidal wave 8,25 . Also, far from the origin ( x = 0 ), it is assumed that there is no concentration exchange with the system. Hence the initial and boundary conditions are: where P 0 and X 0 are the initial uniform concentrations for pollution and DO as x → ∞,P 1 and X 1 are the uniform pollution and DO concentrations respectively when η 1 = η 2 = P 0 = X 0 = 0 , P 2 and X 2 are the average added pollution and DO concentrations respectively, P 3 and X 3 are the amplitudes of pollution and DO variation respectively, η 1 and η 2 are the initial pollution and DO decay lengths respectively, ω is the frequency of the concentration variations 8 . It is known that the behavior of variation of the pollutant concentration is opposite to that of the DO. Hence, it is convenient to choose the sine function in Eq. (8) and the cosine function in Eq. (11).

Problem description
It is noticed that the term µ * (1 − exp(− * z * )) in Eq. (14) and the term exp[−η 1 * z * ] in Eq. (16) have opposite effects on the pollutant concentration so we only consider the two cases where only one of them is present in the model.
Case study one ( P * 1 = X * 1 = 0). In this case, it is assumed that pollutants are injected along the river at a rate µ * and that its value in the upstream is lower than its value in the downstream which is modelled as Case study two ( * = 0 or µ * = 0). In this case, it is assumed that no pollutants are injected along the river and that the pollutant and DO concentrations are initially an exponentially decreasing function of the space variable (Eqs. 16 and 19) 9 .

Analytical solution for the steady state case
It is well known that analytical solutions are important in validating numerical solutions. However, the exact analytical solution of Eqs. (14 and 21) is very difficult if not impossible especially for the nonlinear case. Also, previous studies indicates that solution of the unsteady case is very tedious 4,6 . To simplify the equations, we will consider the steady state and k * = 0 , Pe P * = Pe X * = Pe and U P * = U X * = U * . Hence, Eqs. (14-21) take the form: The solution of Eqs. (22)(23)(24)(25)(26)(27) may be written as: where www.nature.com/scientificreports/ . We have checked that Eqs. (28 and 29) satisfy Eqs. (22)(23)(24)(25)(26)(27).

Numerical solution
Analytical solutions of ADE with limited initial and boundary conditions have very few applications and are very lengthy 26 . The added difficulty of the non-linearity makes obtaining exact analytical solutions extremely hard, if not impossible. Numerical methods do not have such limitations, especially for arbitrary conditions. Thus, numerical solutions are obtained using a two-level explicit finite difference scheme. The truncation error is O(�T * , (�z * ) 2 ) and can be reduced until the accuracy achieved is within error tolerance by taking suitable values of T * and z * 27 . Hence, step-sizes Δ z * =1 and Δ T * =0.002 along the z * -domain and T * -domain, are chosen respectively. The explicit finite difference method is employed to solve Eqs. (14)(15) with initial and boundary conditions (16)(17)(18)(19)(20)(21). The central difference scheme was used for ∂ 2 P * ∂z * 2 , ∂P * ∂z * , ∂ 2 X * ∂z * 2 and ∂P * ∂z * and a forward difference scheme for ∂P * ∂T * and ∂X * ∂T * , with these substitutions, Eqs. (14)(15)(16)(17)(18)(19)(20)(21) can be written as: where i and j refer to the discrete step lengths Δ z * and Δ T * for the coordinates z * and T * respectively, and z * is the grid dimension in the z * direction and z * ∞ is the distance measured from the origin at which ∂P * ∂z * → 0 and ∂X * ∂z * → 0. Equations (30 and 31) represents a formula for P * i,j+1 and X * i,j+1 at the (i, j + 1) th mesh point in terms of known values along the jth time row. www.nature.com/scientificreports/  www.nature.com/scientificreports/

Results and discussions
The numerical solutions of Eqs. (30-37) are illustrated in Figs. 1, 2 and 3 for the common input data: K * 1 = 1 ; K * 2 = 0.5 ; S * = α 1 * = 2 ; η 1 * = η 2 * = 0.98 ; ω * = 2π ; * = 0.6 ; R 0 = 1 ; P * 0 = 0.1 ; P 1 * = 0.5 ; P * 3 = 0.5 ; X * 0 = 6 ; X 1 * = 2 ; X * 3 = 0.5 ; 0 ≤ z * ≤ 70. Figures 1 and 2 show the variations of P * and X * with T * for T * = 0.6, 1.6, 2.6 for case studies two and one respectively for the values U * = Pe P * = Pe X * = k * = 1 . It is clear that as T * increases the value of P * increases and the value of X * decreases. This result agrees with that obtained by 6,13 . These figures emphasize the fact that the effect of the term µ * (1 − exp(− * z * )) is apparent in making the corresponding values of P * in Fig. 2 bigger than in Fig. 1 especially far from the origin 10 . One application of this model is that it predicts the pollutant concentration in the future, which is essential in predicting the following two scenarios: (a) To determine with accuracy the intervals along the river with the highest pollution to focus on its remediation by effective methods at the suitable time. (b) To determine with accuracy the intervals along the river with the lowest pollution from which to withdraw clean water for different uses.
Since k * represents the half-saturated oxygen demand concentration for pollution decay, a direct consequence is its dominant effect on P * while its effect on X * is negligible. Figure 3 shows the variations of P * and X * with k * for k * = 0, 0.3, 1 for case study one for the values T * = 1.6 ; U * = Pe P * = Pe X * = 1 . It is obvious that as k * increases, the value of P * also increases and the value of X * increases slightly but is almost unchanged. This result agrees with that obtained by 7,13 . Figure 4 shows a comparison between the steady state solution given by Eqs. (28,29) and the steady state numerical solution of Eqs. (30, 37) (i.e. when the change with time T * is negligible, which is at T * ∼ = 60 for the following parameter values) and an unsteady state numerical solution (for T * = 20 ), for the parameter values: A very good agreement is found between the steady-state solutions. It is seen that as T * increases in the numerical solution of the unsteady case, the numerical solution approaches the analytical solution in both pollution and aeration. This implies that the numerical solution can be used in both unsteady and steady cases. Thus, the tedious complications of the analytical solution are avoidable.
Since the flux of water Q = A U * (the volume of water crossing A per unit time), the flux Q increases as U * increases which directly leads to an increase in the value of P * while the value of X * decreases. Figure 5 shows the variations of P * and X * with U * for U * = 4, 9, 15 for the steady state solution given by Eqs. (28,29) for the parameter values: Pe = 6 ; K * 1 = 3 ; * = 0.6 ; µ * = 0.5 ; K * 2 = 0.5 ; S * = α 1 * = 2 ; 0 ≤ z * ≤ 70 . It is clear that near the source ( z * = 0 ), as U * increases the value of P * increases and the value of X * decreases. This result agrees with that obtained by 11,13 . www.nature.com/scientificreports/   www.nature.com/scientificreports/ It is known that the Peclet number represents the ratio between the advective transport rate and the diffusive transport rate. Consequently, while the value of the Peclet number remains less than one, the increase in the Peclet number leads to a decrease in the effect of diffusion and thus a decrease in pollutant concentration values. Figure 6 shows the variations of P * and X * with Peclet number for Pe = 0.01, 0.1, 1 in the steady state solution given by Eqs. (28,29) for the parameter values U * = 4 ; K * 1 = 3 ; * = 0.6 ; µ * = 0.5 ; K * 2 = 0.5 ; S * = α 1 * = 2 ; 0 ≤ z * ≤ 40 . It is seen that, near the source ( z * = 0 ), as the value of the Peclet number increases, the value of P * decreases and the value of X * increases. This result agrees with that obtained by 13,28 . Figure 7 shows the variations of X * with µ * for µ * = 12, 14, 16 for the steady state solution given by Eqs. (28 and 29) for the parameter values: Pe = 6 ; K * 1 = 3 ; * = 0.6 ; U * = 9 ; K * 2 = 0.5 ; S * = α 1 * = 2 ; 0 ≤ z * ≤ 70 . As expected, as µ * increases the values of X * decrease along the river. This result agrees with that obtained by 4,7 .
Since the fish will survive when X * > 30%S * , it is obvious that this restricts the value of µ * for a fixed interval of the river at a given time T * . From Fig. 7, it is clear that for the given parameters, the region at which fish survive depends on the value of µ * as follows: for µ * = 16, the fish can survive only in the region z * ≤ 20 , for µ * = 14, the fish can survive only in the region z * ≤ 35 and for µ * = 12, the fish can survive at any region along the river. Thus, the suitable interval along the river and the time for fishing can be determined with high accuracy. The importance of this figure is that it graphically explains the actual procedure, which is directly related to real-life situations requiring predictive tools to help with decision-making. The steps for determining the ideal place and time for placing fish farms are suggested as follows: (1) Taking measurements and calculating all the constant parameters of the model, such as ( Pe P * , S * , α 1 * ,…). The value of µ * obtained in this way is the maximum pollutant concentration that is allowed to be injected along the river for fish survival at a given time. Furthermore, for a fixed interval along the river, we see the change in X * (DO concentration) with respect to the different values of µ * which is essential in determining the ideal place for placing fish farms. It is seen that the procedure and the model are neither economically expensive nor have any environmental risk, and due to their simplicity, they could be put into practical experiments. Table 1 presents the analytical solution in the steady-state case taken as a reference and shows a comparison between different step sizes in the finite difference method. We see from the table that 0.002 is the most suitable value for Δ T * . Both Δ z * = 0.5 and Δ z * = 1 give acceptable results relative to the analytical solution.
Furthermore, Table 1 and Fig. 8 illustrate a comparison between the steady-state analytical solution and the finite difference method and NDSolve via Mathematica. From the table, it is clear that the finite difference scheme is more accurate than NDSolve in this case. In conclusion, a good agreement was found between the three solutions, which is clear in Fig. 8.
It should be noted that when the amplitudes of sine or cosine functions (i.e., P * 3 or X * 3 ) are large in relation to the other parameters of the model, it results in a negative value of the concentrations, which deviates from the real-life meaning.
A final comment on the comparison between case studies one and two: While case study two generalizes previous models, case study one presents a preferable model to get useful information due to the presence of µ * representing the amount of injected pollutants along the river. Thus, obtaining the maximum threshold which acts as guidance towards making necessary decisions.   www.nature.com/scientificreports/

Conclusion
• A generalized pollution-aeration model was established by using time as well as spatial-dependent parameters. Also, general boundary constraints were considered in the two case studies, and special cases of the proposed model agree with its predecessors, making it advantageous to be used in global cases. • The developed model was solved analytically in the steady state case and numerically in the unsteady state case using the finite difference method. A comparison is made between the analytical solution in the steady state case and the numerical solution in its steady state. Also, Mathematica software was used via NDSolve for a comparison with the finite difference method as well as the analytical solution in the steady-state case.
In conclusion, a very good agreement is found. • Case study one, i.e., when considering the injected pollutants along the river, is beneficial in decision-making regarding the regulation of industrial and agricultural practices and, if needed, restricting them. Since tons of fish die every year in the Rosetta Branch of the Nile River in Egypt and other rivers suffer from the same problem as well, this study gives a suitable solution to this problem. The suitable time for fishing as well as the location of fish farms are determined with high accuracy by following a simple step-by-step procedure.

Limitations and recommendations
This work is limited by the usage of a conditionally-stable numerical method, namely, the finite difference method. The analysis carried out revealed that for the given parameter values, the results are satisfactory. However, using an unconditionally-stable numerical method such as the Crank-Nicolson method for example, will allow for even wider parameter values. Hence, analyzing the general aeration model using different numerical methods could be a fruitful area of research. Also, considering different boundary conditions such as pulse-type and mixed-type conditions is recommended.